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Abstract. Motivated by the experimental ability to produce monodisperse particles in microfluidic devices, 
we study theoretically the hydrodynamic stability of driven and active crystals. We first recall the theo- 
retical tools allowing to quantify the dynamics of elongated particles in a confined fluid. In this regime 
hydrodynamic interactions between particles arise from a superposition of potential dipolar singularities. 
We exploit this feature to derive the equations of motion for the particle positions and orientations. After 
showing that all five planar Bravais lattices are stationary solutions of the equations of motion, we consider 
separately the case where the particles are passively driven by an external force, and the situation where 
they are self-propelling. We first demonstrate that phonon modes propagate in driven crystals, which are 
always marginally stable. The spatial structure of the eigenmodes depend solely on the symmetries of the 
lattices, and on the orientation of the driving force. For active crystals, the stability of the particle positions 
and orientations depends not only on the symmetry of the crystals but also on the perturbation wave- 
lengths and on the crystal density. Unlike unconfined fluids, the stability of active crystals is independent 
of the nature of the propulsion mechanism at the single particle level. The square and rectangular lattices 
are found to be linearly unstable at short wavelengths provided the volume fraction of the crystals is high 
enough. Differently, hexagonal, oblique, and face-centered crystals are always unstable. Our work provides 
a theoretical basis for future experimental work on flowing microfluidic crystals. 



1 Introduction 

The dynamics of passive suspensions is a field with a long 
history in physical hydrodynamics. Much effort has been 
devoted to understand, e.g., the origin of fluctuations in 
the sedimentation of spheres under gravity as well as in- 
stabilities in suspensions of elongated fibers (see reviews 
in |1I2| and references therein). More recently, a signifi- 
cant experimental |3I4I5| and theoretical [61718] research 
effort has focused on the dynamics of active suspensions 
where instead of having particles driven by an external 
field (e.g. gravity), one considers the dynamics and inter- 
actions of self-propelled synthetic or biological swimmers. 
In this case, the interplay of activity and hydrodynamic 
interactions leads to long- wavelength instabilities [9110] . 

Most of the past work on passive (driven) and ac- 
tive suspensions has focused on instabilities and fluctuat- 
ing behavior in three-dimensional systems. However, over 
the last ten years microfluidics has offered a number of 
simple and effective solutions to produce and manipu- 
late large ensemble of highly monodisperse microparticles, 
prone to form crystal structure in quasi-two dimensional 
channels [11]. For driven particles, these technological ad- 
vances have motivated, for example, the study of the non- 
linear dynamics of finite flowing crystals |12ll3j , phonons 
in one dimensional microfluic-droplet crystals [14 and flow- 



ing lattices of bubbles [TS]. In the case of active parti- 
cles, these fabrication methods could be extended to self- 
propelled catalytic colloids |16|17j or reactive droplets [18] . 



Motivated by these advances, we take in this paper an 
approach contrasting with the traditional study of disor- 
dered suspensions and consider the dynamics of confined 
driven and active hydrodynamic crystals. We first develop 
a formalism to study theoretically position and orienta- 
tion instabilities for flowing discrete suspensions under 
confinement. In the case of driven particles, we demon- 
strate formally that all crystals are marginally stable and 
study in detail the eigenmodes of deformation for all five 
two-dimensional Bravais lattices. For active particles, we 
show that square and rectangular crystals are linearly un- 
stable at short wavelengths provided the volume fraction 
of the crystals is high enough. Differently, hexagonal and 
oblique (resp. face-centered) crystals are always unstable 
for long- (resp. short-) wavelength perturbations. In con- 
trast with past work on three-dimensional swimmers sus- 
pensions, the stability of confined active crystals is found 
to be independent of the pusher vs. puller nature of the 
actuation of individual active particles |19) . 
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Fig. 1. Schematic representation of the problem ad- 
dressed in this paper: An extended hydrodynamic crystal 
is composed of anisotropic particles confined in a channel 
of height h which are either actively swimming or passively 
driven by an external force (top- and side-views). 

2 Theoretical setup 

2.1 Particle crystal in a confined fluid 

We start by describing the theoretical framework we use 
to quantify the large scale dynamics of both active and 
driven microfluidic crystals. We focus our study on the 
case of identical particles living in quasi-bidimensional flu- 
ids, as sketched in Fig. [T] The fluid is Newtonian and has 
a homogeneous thickness h in the z-direction, comparable 
to the size of the particles. Our formalism will be valid 
both for thin films lying on a solid substrate (with one 
free surface and one no-slip wall), as well as for microflu- 
idic geometries where the fluid is confined between two 
parallel plates. The particles can be either axisymmetric 
or anisotropic and are organized in two-dimensional crys- 
tal, see Fig. [I] If a denotes the typical lattice spacing of 
the crystal and b the typical extent of the particle in the 
(x,y) plane, we consider in this paper the dynamics in 
the dilute limit, e.g. a 3> b. In this limit, each particle i is 
appropriately modeled as a pointwise body characterized 
by its in-plane position, Rj(i) = (xi(t),yi(t)), and its in- 
plane orientation, Pi(i), where pi is a unit vector making 
an angle 9i(t) with the x-axis. Having microscopic systems 
in mind, we neglect the particle inertia and work in the 
limit of zero Reynolds number. In this Hele-Shaw setup, it 
is a classical result that the fluid flow is potential |20] . The 
2-averaged fluid velocity, V and the z-averaged pressure, 
P, are therefore related by 

V(r) = -GVP, (1) 

where G = ah 2 /rj; here rj is the fluid viscosity, and a = 1/3 
for a thin film, and a = 1/12 for a shallow microchannel. 
Together with incompressibility, V • V = 0, Eq. [I] deter- 
mines the fluid flow and stress away from the particles. 



We henceforth consider either swimmers moving along 
their principal axis p^, or passive particles driven by a 
uniform force field oriented along the x-direction (gravi- 
tational, electrostatic, magnetic,...). In all cases, the speed 
of an isolated particle in a quiescent fluid is constant and 
denoted Uq . In addition to their individual dynamics, par- 
ticles also follow the surrounding flow, and the equation 
of motion for particle i thus reads 

a t R,-c/ q + /iV(R i ), (2) 

where q = Pi for swimmers, q = x for driven particles, 
and fi is a non-dimensional mobility coefficient [14121) . 
Passive tracers have \i = 1. Conversely, for thick particles, 
the friction against the solid wall(s) can significantly re- 
duce the advection speed, which is smaller that the local 
fluid velocity, and thus < fx < 1. In principle, /i should be 
a tensor for anisotropic particles but for simplicity we con- 
sider only particles which are weakly anisotropic and thus 
/i is assumed to remain a scalar Q In addition to a change 
in their velocity, anisotropic particles experience hydrody- 
namic torques which favor an orientation along the local 
elongation axis of the flow. This classical hydrodynamic 
result, which can also be anticipated from symmetry ar- 
guments, leads to the so-called Jeffery's orbits 22J. As 
the flow is irrotational (potential flow), the orientational 
dynamics reduces to 

dtpi = 7(J-p i pi)-E(R i )-p i , (3) 

where E is the strain rate tensor, E = ~ [VV + (VV) T ] , 
and 7 > is a rotational mobility coefficient which is non- 
zero for anisotropic particles and zero for axisymmetric 
bodies. 



2.2 Long-range hydrodynamic interactions 

As a particle located at R(t) moves in the fluid, it in- 
duces a far-field velocity, denoted v(r — R), at position r. 
A given particle i responds to the flow induced by all the 
other particles in the crystal, and is therefore advected at 
velocity V(Rj) = MX)«i v (R-i ~ ^-j)- I n this section we 
provide a quantitative description of the far-field hydro- 
dynamic coupling, v, between the particles; for the sake of 
clarity, we separate the case of passive and active particles. 

2.2.1 Hydrodynamic interactions between driven particles 

In the driven case, each particle of the crystal is subject to 
a constant external force, f = /x, which results in a far- 
field perturbation which we denote v 1 and is the Green's 
function of Eq.[T] In three-dimensional flows, the response 
to a force monopole is known as a Stokeslet, and de- 
cays spatially as ~ 1/r. In our quasi-2D geometries, solid 
walls act as momentum sinks and screen algebraically the 

1 Note that the anisotropy of the mobility coefficient is much 
weaker in confined than in unbounded fluids due to the short 
range of hydrodynamic interactions in quasi 2D geometries. 
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Fig. 2. Sketch of the dipolar flow field (potential source 
dipole) induced by driven particles (left) and active swim- 
mers (right). 

Stokeslet contribution, which then decays as v 1 ~ 1/r 2 
and takes the functional form of a potential source dipole, 
as shown in [23 23]. In addition, the particles have a fi- 
nite size and their advection by the surrounding fluid is 
hindered by the lubrication forces induced by the confin- 
ing walls (even in the absence of external driving). Due 
to incompressibility, any relative motion with respect to 
the fluid results in another algebraic far field contribution, 
v 2 . As shown, e.g. in [14], v 2 has also the form of a poten- 
tial dipole with the same spatial decay, v 2 ~ 1/r 2 . (We 
note that in unbounded fluids, this potential contribution 
scales as 1 /r 3 and is thus subdominant with respect to the 
flow induced by a pointwise force, which decays as 1/r). 
Therefore, in confined flows, the two contributions have 
the same form |14|23j and the overall flow disturbance, 
v d = v 1 + v 2 , takes the form of a z-dipole, 



>(r) 



27rr 2 



(2rr - I) ■ x, 



(4) 



where r = |r| and the dipole strength, cr, is the sum of 
the two contributions, a — Ab 2 Gf + Bb 2 Uo, where A and 
B are two dimensionless shape factors (/ is the identity 
tensor) . The symmetry of the streamlines for this flow field 
are illustrated in Fig. [5] (left). 

2.2.2 Hydrodynamic interactions between active swimmers 

By definition swimmers do not require an external force to 
propel themselves. The stress distribution on the surface 
of a self-propelled particle has thus, at least, the sym- 
metry of a force dipole [25]. The canonical theoretical 
setup used to describe (dilute) suspensions of swimmer is 
to consider an ensemble of such force-dipoles as all other 
multipolar contributions to the far field are subdominant 
in an unbounded fluid [9 19 26 . However, as mentioned 
above, confinement results in an algebraic screening of the 
hydrodynamic interactions. In the quasi-2D geometry at 
the center of our paper a force dipole decays spatially as 
~ 1/r 3 , a contribution which is therefore subdominant 
compared to the ~ 1/r 2 potential dipole arising from in- 
compressibility (similarly to driven particles) [53] • For ac- 
tive swimmers, the far-field flow disturbance has thus also 



the symmetry of a potential source dipole, the difference 
with the passive case being that the dipole direction is 
now the swimmer orientation (Fig. [2j right). For a swim- 
mer orientated along p, we obtain a flow given by 



v s (r,p) = 



27rr- 



(2ff - I) ■ p, 



(5) 



with the dipole strength a — Bb 2 Uo (B is the same shape 
factor as in Eq. [4]). We therefore see that, that in con- 
fined fluids, the usual distinction between pushers and 
pullers swimmers (contractile and extensile), which is at 
the heart of qualitatively different behaviors in unconfined 
fluids |7I8| , is irrelevant. The magnitude and sign of the 
induced dipolar flow is solely set by that of the swimming 
speed, irrespective of the microscopic swimming mecha- 
nism. 

In summary, Eqs. p] pi and either Eq. H] (in the driven 
case) or [5] (active case) fully prescribe the dynamics of 
the discrete particle positions and orientations. As noted 
above, the main difference between active and passive par- 
ticles concerns the orientation of the dipolar flow field: the 
orientation is slaved to the swimmer direction for active 
particles whereas it is constant and aligned along the x- 
dircction for driven particles (the difference is further illus- 
trated in Fig. [2j. We will show in the next sections that 
this distinction markedly impacts the large-scale crystal 
dynamics. 



3 Are hydrodynamic crystals stationary? 

When addressing the dynamics of an ordered phase, the 
first important question is whether this phase does corre- 
spond to a stationary state. We focus here on the five pla- 
nar Bravais crystals, which encompass all possible symme- 
tries for bidimensional mono-atomic crystals (see Fig. 3j. 

Let us first consider the case of crystals composed of 
driven axisymmetric particles. The equations of motion 
reduce to 



c\R 4 = (7 x + A t^v d (R i 



R 7 " 



(6) 



where the Rj's belong to one of the Bravais crystals from 
Fig. [3] The lattice structure is conserved provided that 
dt(Ri — Rj) =0 for all i and j. It follows from Eq. [6] that 

a t (R t -R J )= A1 E fe ^v d (R I -R fc )- At E^ 5 v d (R J -Rfe)- 



By definition, all crystals are invariant upon translation 
along (Rj — Rj) , which readily implies that the two sums 
are equal, and therefore that any driven crystal made of 
axisymmetric particles is stationary structure (i.e. dt(Ri — 
R,-)=0). 

To extend this result to driven crystals composed of 
anisotropic particles we first recall that all the Bravais lat- 
tices are invariant upon the parity transformation r — > — r. 
Moreover, as v d (r) is invariant upon this transformation 
whereas the sign of the gradient operator is reversed, we 
see that the strain rate tensor constructed from a superpo- 
sition of potential source dipoles, Eq.[4] transforms accord- 
ing to E(— r) = — E(r). This implies that for any Bravais 
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Fig. 3. Geometry of the five planar Bravais lattices. Anisotropic cells are characterized by the ratio, e, between the 
two cell dimensions. The angle j3 is the tilt angle of the oblique and hexagonal cells. For each particle labeled "0" we 
also display and number all nearest neighbors. 



crystal, E(Rj) is identically zero anywhere on the lattice. 
It follows from the equation for the orientational dynam- 
ics, Eq.[3j that dtPi = for all the particles. In conclusion, 
in the driven case, both the crystalline structure and the 
particle orientations remain stationary (in other words, 
the crystals are fixed points of the dynamical system). 

It is straightforward to generalize the above results to 
swimmer crystals. The equation of motion for the posi- 
tions, Eq. [6j is given by 



d t R t = U pi 



Rj , p j / 



(7) 



Obviously, d t (ELi — Rj) cannot be zero if the f>i are not 
all identical. Therefore, the crystal structure cannot be 
conserved if the initial orientation of the particles is not 
uniform - in such cases the crystal would "melt" . For uni- 
form orientations, say along x, Eqs. [6] and [7] are identical, 
and so is the equation for the orientational dynamics since 
v s (r, pi) — v d (r). We are thus left with the same problem 
as in the driven case, which implies that the structure of 
the crystals is conserved as long as the particles all swim 
along the same direction. 



4 Driven hydrodynamic crystals are 
marginally stable 

We start by investigating in this section the linear stability 
of the five Bravais crystals with respect to perturbations in 
both the position and the orientation of the particles, with 
a special focus on the experimentally-relevant square and 
hexagonal lattices. Anticipating on our results, we note 
that the geometrical classification in term of the Bravais 
lattices might not necessarily be relevant to the dynamics 
of flowing crystals. 

In order to proceed we make use of two additional as- 
sumptions. Firstly, we consider the case of particles uni- 
formly aligned along the x-axis prior to the perturbations, 
as depicted in Fig. [3j Secondly we assume that the driv- 
ing force is aligned with one of the principal axes of the 
crystal. Our following study can be easily extended to a 
more general setup, but this would make the formula and 
the discussions much more tedious. 

We denote JR.; and Spi ~ (^y the infinitesimal pertur- 
bations of the particle positions and orientations respec- 
tively, so that Ri — > Hi + <5Rj, and — > x + 6$. Using 



the property that E = for dipoles organized into a Bra- 
vais lattice (as discussed in the previous section) , and after 
some algebra, the linearization of the equations of motion, 
Eqs. |3]and[6j yields 



d t SR, = [Vv d (Ry)] • SRij, 



and, 



<V4(Ry)]HRi. 



(8) 



(9) 



where Rjj = Ri — Rj, and SRij = 6R4 — SRj. Eqs. [8] 



and [9] dictate the dynamics of the elementary excitations 
in the frame where the unperturbed crystal is stationary. 
We note that the direction of the crystal translation is, in 
general, different from the driving direction. 

We now exploit the symmetries of the dipolar inter- 
actions. Inspecting the flow given by Eq. [4j we deduce 
that d x v d = dyV d , and d x v x = —d y Vy. Using these rela- 
tions, we look for plane waves solutions, (6Xi,6Yi,6i) = 
(8X,SY,0) exp(iujt — iq • R,-). By doing so, we obtain a 
linear-stability system, which we write in the generic form 




Mi Mi 
Mi -Mi 
M 3 M 4 



where the coefficients of the stability matrix M are 




(10) 



En 


— exp(iq 


Ry)] a^Rij), 


(11) 


En 


— exp(iq 


Ry)] 5 x ^(R y ), 


(12) 


En 


— exp(iq 


Rij)] d xx v y (ILij), 


(13) 


En 


— exp(iq 


Rii)] d yx v*(Rij). 


(14) 



We readily deduce from the matrix structure that a 
perturbation in orientations only would not induce any 
change in the crystal conformation. This is a direct conse- 
quence of the dipolar coupling between the particles, v d , 
which is only a function of the driving force direction and 



Nicolas Desreumaux et al.: Active and driven hydrodynamic crystals 



5 



not of the particle orientation. On the contrary, pertur- 
bations in the position of a particle modify both position 
and orientation. In addition, perturbations in orientation 
only neither relax, grow or propagate. As the third col- 
umn of the matrix M is always 0, this implies that it will 
always admit the eigenvalue loq — 0, associated to the 
pure-orientation eigenmode (0,0, 1). 

The other two eigenvalues of the M matrix are uj± = 
± y/ Mf + M| . Exploiting again the fact that all Bravais 
lattices are invariant upon parity transformation to write 

Mi = ]T [1 - exp(zq • Ry)] a^(Ry) 



I f I. 



^[l-exp^q-R^jS^CR^). (15) 



By definition R.y = — Rjj, and due to the dipolar sym- 
metry of the hydrodynamic interaction, we ha ve d x v x (r) = 



15 



we infer 



-d x v x (— r). Using these two equalities in Eq. 
that Mi — — fx Y^j sin(q ■ ILij)d x v x (Tiij), and therefore 
Mi is a real number. Using the same method, and the 
identity d x v y (r) = — d x Vy(— r), one can show that M2 
is real as well. Therefore, for any symmetry of the crys- 
tal, the pulsations of the plane waves, u>±, are real. In 
other words, for any Bravais lattice the crystal structure 
of driven particles is dynamically marginally stable. 

Notably, the dipole strength a can be eliminated from 
the equations of motion Eqs. [5] and [9] by rescaling the 
timescale. Therefore, the linear stability of the monocrys- 
tals is a purely geometrical problem. The corresponding 
eigenmodes do not depend on the translational speed Uq, 
but only on the orientation, and on the symmetries of the 
lattice. 

Interestingly, we see that phonons propagate with the 
pulsations w±, despite the fact that particles have no in- 
ertia and that no potential forces couple the particle dis- 
placements. This seemingly counterintuitive result gener- 
alizes the experimental observations made by Beatus and 
coworkers in [14] where they revealed that sound modes 
propagate along ID-droplet crystals flowing in quasi-2D 
microchannels. These results are, importantly, specific to 
the quasi-2D geometry, which is relevant for numerous mi- 
crofluidic and thin films applications. In unbounded fluids, 
the change in the symmetry of the hydrodynamic interac- 
tions results in the destabilization of the crystal structure 
as shown theoretically and experimentally |27j . 

Below, we derive the dispersion relation for each of the 
five Bravais crystals, with a special attention given to the 
case of square and hexagonal lattices 



4.1 Square lattice 

In order to compute the coefficients of the M-matrix an- 
alytically we now make a nearest neighbor approxima- 
tion. In a similar context, this approximation has proven 
to yield qualitatively correct results for unbounded fluids 
[2"7] . We introduce a reference particle labelled as 0. The 
four nearest neighbors in the square crystal are labeled as 




Fig. 4. Normalized dispersion relation for the square lat- 
tice plotted from Eq. [l7Kw + only), for 2\iaj'Kc? — 1, and 
a=l. 



1, 2, 3, and 4 (Fig 
the coefficient of t 



3]) . In this geometry, we easily compute 
ic M matrix as 



M = 



2/i<7 



7ra° 



sin(q x a) 
- sm(q y a) 




372 
1" 



— sin(q y a) 

- sm(q x a) 
^[cos(g x o) + cos(q y a) 



<r 


2] 
(16) 

The dispersion relation of the infinitesimal excitations 
can be deduced by diagonalizing M. The three eigenvalues 
are uiq — 0, and 



LU± = ± 



2/1(7 



^ \J b\vl 2 (q x a) + sin 2 (q y a). (17) 



This dispersion relation is plotted in Fig. [4] 

To gain insight into the propagating modes, we focus 
on the large-scale (long wavelength) response of the crys- 
tals. Expanding Eq.[l6]at leading order in the wave vector 
amplitude for q —¥ 0, we find 



M = 



2/i er 



q x -q y 
-q y -q x 




0(q 2 ). 



(18) 



We notice that in this small-q limit, the orientation 
and the position degrees of freedom are totally decoupled. 
The three eigenvalues are u>o — 0, and u>± = ±2(iaq/ '(wa 2 ). 
The two non-trivial modes are non-dispersive and propa- 
gate with a constant "sound-velocity" c± — ±2(xa/{j:a 2 ), 
which increases with the magnitude of the hydrodynamic 
coupling. 

To understand physically how the excitations propa- 
gate, we focus on two specific cases. Let us first consider 
longitudinal perturbations, q = qx. The mode w_ is here 
associated with the eigenvector (0, 1,0). It corresponds to 
shear waves which propagate in the direction opposite to 
the driving, as illustrated in Fig. [5|\. The second sound 
mode (w + ) corresponds to compression waves along the 
x-axis propagating in the driving direction, see Fig. [5j3. 
The corresponding eigenvector is (1,0,0). 
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Fig. 5. Sketch of the propagative eigenmodes in a square 
lattice for q y = 0. The full line corresponds to the direction 
of the driving, the dotted line indicates the direction of 
the wave propagation. A: Shear modes; B: Compression 
modes. The particle orientations are not affected by the 
perturbation. 



For excitations propagating in the direction transverse 
to the driving, q = qy, the eigenmodes couple the dis- 
placements along the two principal axes of the crystal. 
The mode w_ is associated with the eigenvector (1, 1, 0). 
It corresponds to the superposition of a compression mode 
in the y direction, in phase with a shear in the x direc- 
tion. The second mode (w+), with eigenvector (—1,1,0), 
is a combination of a dilation in the y direction, which 
propagate in antiphase with a shear wave in the x direc- 
tion. 

To close, we note that the dispersion relation of the 
phonons remains unchanged if the driving force is not 
aligned with one of the principal axes of the crystal, al- 
though in that case the form of the eigenmodes is more 
complex. 



4.2 Hexagonal lattice 

We now consider the case of the hexagonal lattice. The 
main technical difference with the square lattice is that 
the reference particle has now six nearest neighbors, see 
Fig. [3] Repeating the same procedure as above, we com- 
pute the coefficients of the stability matrix and obtain 



M 



2/1(7 

7TQ, 3 



' M[ s 
-M[ 
M3 M4 



(19) 



where 





= sin(g x 




37 y3i 








= -Mr, 

jjLO, L 



tq y v3a \ 
" 2 /' 



(20) 



Notably, the upper-left 2x2 sub-bloc of M is diagonal. 
As a consequence, an excitation of the position along one 



& 




Fig. 6 

lattice plotted from Eq 
and a = 1. 



Normalized disp ersi on relation for the hexagonal 
(lu + only), for 2/ia/ira 2 



21 



1. 



direction (x or y) induces no net displacement in the trans- 
verse direction. The dispersion relations, plotted in Fig.[6j 
is given by 



u>± = ± 



2/1(7 



sin(q x a) — 2 sin( 



q y V3a 
2 ' cos (^-) 



q x a 



(21) 



We see from Eq. 21 that there exist two specific orien- 
tations of the wavevectors for which no excitation prop- 
agates. For perturbations making angles equal to 7r/6, 
(resp. 7r/2) with the x-axis, we obtain M[ = (resp. M[ = 
M3 = 0); in these cases, the matrix is not diagonalizable 
and the only solutions of Eq. 19 is u — 0. Phonons there- 



fore do not propagate in those two directions. 

To illustrate the difference in the dynamics between 
the hexagonal and the square crystals we consider the be- 
havior in the long wavelength limit. Expanding Eq. |19| at 
leading order in q, we obtain 



M = 



2/1(7 







97 • 

4fia 2 l 1x% 



M«2 



2 «S 



(22) 



and the eigenvalues takes the form u)± — ±^\3q x qy — q x \- 
We infer from this formula that hydrodynamic crystals 
having an hexagonal symmetry are "softer" than square 
crystals. When q — > 0, the speed of sound goes to as 
q 2 and even the large- wavelength phonons are dispersive. 
Furthermore, the sound modes couple the displacements 
and the orientation of the particles. 

To convey a more intuitive picture, we again focus on 
two specific directions of propagation. We first consider 
longitudinal perturbations along the first principal axis 
of the crystal, q = gx. As above, we find that one of the 
eigenvector correspond to a pure compression along the x- 
axis, (1, 0, 0). The second eigenvector, (0, —iq/9, 1), mixes 
shear and orientational waves (bending modes) in quadra- 
ture, and depends explicitly on q. Such a coupling was not 
observed for the square lattice. A second simple case con- 
cern the excitations propagating along the second princi- 
pal direction, namely q = cos(7r/3)x + sin(7r/3)y. Here, 
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Fig. 7. Normalized dispersion relations for the rectangu- 
lar, oblique and centered rectangular lattices (w+ only), 
for 2fi(r/TTa 2 = 1, and a = 1. Rectangular: e = 0.8. 
Oblique: e = 0.9 and /3 = (tt/2) - 0.2. Centered rect- 
angular: e — 0.9. 

the two eigenmodes mix the particle displacements (in 
only one of the two directions, since x and y cannot cou- 
ple) and their orientation. They are given by (0, —2iq/9, 1) 
and (— 2zq/(9v / 3), 0, 1) and correspond to w_ and oj + re- 
spectively. 

4.3 Rectangular, oblique, and face-centered lattices 

The lattice geometries corresponding to the rectangular, 
the oblique and the face-centered lattices are shown in 
Fig. [3] To derive the eigenmodes we restrict our analysis 
to calculations with four nearest neighbors, an assumption 
which restrains the number of crystals for which our cal- 
culations are correct (weakly anisotropic and weakly tilted 
lattices only, as sketched in Fig. [3]) . It is straightforward 
to proceed mathematically as in the two previous cases 
and derive the two sound modes, uj±, propagating. The 
results for the dispersion relation are plotted in Fig. [7j 
In the small g-limit, these phonons always propagate in a 
dispersive manner. As q goes to zero, the sound velocities 
reach a constant value which depends on the orientation 
of the propagation due to the crystal anisotropy. 

4.4 Response of driven hydrodynamic crystals to finite 
amplitude perturbations 

Before closing this section, we make a final remark re- 
garding the stability of all the five Bravais crystals with 
respect to finite amplitude perturbations. We start by a 
simple observation on the relationship between the vari- 
ous lattices. A rectangular lattice corresponds to a square 
lattice transformed upon a finite homogeneous stretching. 
An oblique lattice is obtained by stretching and shear- 
ing a square lattice. The hexagonal lattice is an oblique 
lattice with a tilt angle of 7r/3. Finally, a face-centered 
lattice is obtained from a rectangular lattice by apply- 
ing a shear modulated at the highest possible wavelength 
(q = 2-k/o). As all these structures are stationary, and 
marginally stables at the linear level, we can deduce that 
any finite amplitude deformation corresponding to an ho- 
mogeneous shear, or stretch, of the crystal would also be 
a marginal perturbation: their growth rate would be zero. 
The same conclusion also holds for rectangular crystals de- 
formed by the specific high-g shear that would transform 
them into a face-centered lattice. 



5 Hydrodynamic stability of active crystals 

We now move on to investigate the linear stability of active 
swimmer crystals. To do so, we use the same theoretical 
framework as in the previous section. The swimmers self- 
propel along one of the principal axes of the crystals. We 
also recall that in this active case, the swimming direction 
is slaved to the particle shape and so is the dipolar flow 
(Fig. [2]). Following the same strategy as in the case of 
driven particles, we first establish the linearized equations 
of motion. Combining Eqs. [3j [5] and [7] we obtain 

dtSR, = U e t e y (23) 
+ fiJ2 ([Vv s (%,x)] ■ SRii + [9 e v s (R 4J ,x)] 6-) , 

and 

dA = 7 ]T [V [d x v s y (Ry , x)] • SRti (24) 

+ d e E s yx (R ij> p j =x)e j }, 

where E* is the (y, x) component of the strain rate tensor 
associated with the dipolar perturbation induced by the 
swimmer located at j, namely v s (Rjj,pj = x). 

Two important remarks can be made at this point. 
First we see that the stability equations now depend ex- 
plicitly on the swimming speed of the particles. In addi- 
tion, contrary to driven lattices, the stability of the swim- 
mer crystals depends on the particle shape through 7. 
Therefore, we discuss below isotropic and anisotropic swim- 
mers separately. 

5.1 Isotropic swimmers 

Isotropic particles correspond to 7 = 0. Their dynamic 
equation are significantly simplified as Eq.[24]is now trivial 
and the swimmer orientation remains constant. Since the 
flow is irrotational, no hydrodynamic torque (from vortic- 
ity) is present to modify the orientations of the particles. 
As in the previous section on driven suspensions, we look 
for plane waves solutions, from which we infer the form of 
the stability matrix M. This matrix is analogous to the 
one defined in Eq. [10] but takes here a slightly different 
structure 

(Mi M 2 M 5 \ 
M = M 2 ~Mi M 6 , (25) 
\ J 

where Mi and Mi are given by Eq. fTTJ and Eq. [12] respec- 
tively. The two new coefficients are 

M 5 = [1 - exp(iq • Ry)] fyt£(Ry,x), (26) 

M 6 = -1U0 - ifiJ2 t 1 ~ ex P( ic l ' R 'j)] 9 K(Rij,x).(27) 

Independently of the crystal symmetry, we see that 
the eigenvalues of the above matrix are identical to the 
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Fig. 8. Unstable position/orientation mode for a square 
lattice of active particles if the volume fraction is high 
enough (p < 1). The mode is a compression along the y- 
direction out of phase with a splay perturbation of the 
particles orientation and leads to the formation of short 
wavelength bands. 



lattice. The parameter p quantifies the relative magnitude 
of the swimming speed and the dipolar advection velocity 
induced by a neighboring particle. Recall that a is itself a 
function of Uq, and of the particle shape, and a = BUob 2 , 
where B is a shape factor of order 1. Therefore, p scales as 
p ~ fi~ 1 (a/b) 2 . Large values of p correspond to the dilute 
limit, a 3> b, in which our far field approach is expected 
to be quantitatively correct. Small values of p correspond 
to a dense crystal, for which our model should capture the 
essential physical features. The presence of p in the matrix 
M means that the crystal stability now strongly depends 
on the particle volume fraction. 



one we found for driven crystals, loq — 0, and uj± = 
± \J M 2 + M| . Therefore, crystals composed of isotropic 
swimmers are marginally stable and phonons propagate 
with the same dispersion relations as in driven lattices, 
albeit with different eigenmodes. 



5.2 Anisotropic swimmers 

We now explore the richer phenomenology arising from 
swimmer anisotropy. Generic results cannot be established 
in a framework as general as in the isotropic case. We 
proceed with the calculation under the nearest neighbors 
approximation, and deal with the five Bravais lattices sep- 
arately. 



5.2.1 Square lattice 

To establish the linear stability of the square crystal we 
compute all the coefficients of the M matrix using Eqs.[23| 
and EH and obtain 



2ua { 1 2 
M = M' 2 -M[ Mg , . 

wa \ M' 4 M£ 



(28) 



with 



M[ = sin{q x a), (29) 
M' 2 = s\n{q y a), (30) 

M 4 = — 3i — [cos(q x a) + cos(q y a) — 2] , (31) 

M' & = -ia (p+^ [cos(q y a) - cos^a)]^ , (32) 

(33) 



Mj = - sin(q x a), 

where we introduced the dimensionless number 

irU a 2 



P 



2fia 



(34) 



Unlike the driven case, the stability matrix for active 
particles is not characterized solely by the geometry of the 



As even in the large-p limit the eigenvalues of M take 
a quite complex form, we proceed to consider the short 
and long wavelength excitations separately. In the limit 
q — >• 0, the matrix M has again three real eigenvalues, cor- 
responding to three propagating modes with frequencies 
luq = 2jaq x /(ira 2 ), and U)± — ±2/icrq/(7ra 2 ). The mode 
uq is a combination of phonons and orientation waves, 
whereas lu± are the phonon modes we found for driven 
crystals (Fig. [5]). 



In the high-g limit (small wavelengths) the phenomenol- 
ogy is markedly different. For waves vectors of the edge of 
the Brillouin zone, q x = and q y — it /a, we find uq — as 
well as two non-trivial modes, uj± = ±^^--\/67/i(p — 1). 
Importantly, uj± are either real or pure imaginary num- 
bers depending on the magnitude of p. In principle, p > 1 
for dilute crystals, and therefore the modes uj± corre- 
spond again to a combination of phonons and orienta- 
tion waves. However, we can expect our results to hold 
at a qualitative level for more concentrated systems, for 
which p < 1. In such a case, the hydrodynamic cou- 
pling destroys the square crystal structure. Specifically, 
the w_ mode is unstable. It correspond to the eigenvector 
(5X, SY, 8) = (0, — i^r^u!-, 1), which combines a compres- 
sion along the y-axis and splay distortions of the particle 
orientation. In this strong hydrodynamic coupling limit, 
the square crystal evolves to form short wavelength bands 
aligned with the average swimming direction, as sketched 
in Fig. [8] 



At second order in q — > and given that p is small 
enough, the eigenvalues u>± have a non-zero imaginary 
part which scales as q 2 . These eigenvalues correspond to 
the roots of a 3rd-order polynomial, which has no analyt- 
ical solution. Therefore, we proceed to a numerical inves- 
tigation of the short wavelength dynamics of the crystal. 
We compute numerically the eigenvalues of the matrix M 
for all q's and < p < 3. We find that the square crys- 
tals are indeed always unstable for p < 1. In addition the 
wavenumbers q x — and q y = it / a correspond to the most 
unstable mode as shown in Fig. [9] for p= 1/2. 
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Fig. 9. Normalized growth rate (iuj-) of the (q x ,q y ) 
modes of the square lattice of active particles for p = 1/2. 
The parameters are 2[Mr/ira 2 = 1, a = 1 and j/fi = 1. 
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Fig. 10. Normalized growth rate (ioJ-) of the (q x ,q y ) 
modes of the hexagonal lattice of active particles for p = 5. 
The parameters are 2/Mj/ira 2 = 1, a = 1 and j/fj, = 1. 



5.2.2 Hexagonal lattice 

When the lattice has hexagonal symmetry, the structure 
of M is somewhat simplified and we obtain in this case 



M 



2[ia 

TXO? 



( M[ 



-M[ -iap - ^M'i 
M' 3 Mi IM[ 



with, 



M[ = sva{q x a) — 2 sin(g x a/2) cos(v3g y a/2), 
Ms = i jT^l sm(q x a/2) sm(V3q y a/2), 



.1/; :'./ 



//.a 



(35) 

(36) 
(37) 



cos 



(q x a/2) cos (V3q y a/ 2) — cos(q x a) . (38) 



Similarly to the square lattice, M does depend on the 
relative magnitude of the hydrodynamic coupling through 
p. The general form of the eigenvalues is too complex to 
yield an intuitive picture. However, the salient features 
correspond to small wave vectors. In this limit q ~ > 0, the 

eigenvalues of M are lu q = 0, and u± — \^np{q 2 — q 

For all values of p, there exists therefore an infinite num- 
ber of unstable modes growing at a rate \w±\. They cor- 
respond to perturbations in the position of the particle 
propagating along a direction making an angle comprised 
between [±7r/4, ±37r/4] with respect to x. This behavior is 
illustrated in Fig. [lOj where we show the variations of the 
imaginary part of w± in the (q x ,q y ) plane for p — 5. We 
note that the most unstable mode is again a combination 
of compression along the y-axis and splay-like instability 
of the particle orientation. 



5.2.3 Rectangular crystals 

The behavior for the rectangular lattices is very similar to 
what we found for square lattice (within the limits of the 



nearest neighbor approximation). These crystals are all 
stable for long wavelengths but can display short wave- 
length instabilities. Denoting e the aspect ratio of the lat- 
tice cell (Fig. [3| a numerical diagonalization of the stabil- 
ity matrix reveals that again the most unstable mode lies 
on the edge of the Brillouin zone in the y direction. The as- 
sociated eigenvalue is w± = 7r (^ 1 )3 \/3m7 [(?P ~ l) g2 ~ 
Hence, there exists a critical value p c = |(1 + e~ 2 ) such 
that the crystal destabilizes for p < p c ; note that p c is 
a decreasing function of the aspect ratio which plateaus 
at p = 1/2. Dilute crystals corresponding to high p- values 
are stable and display phonons modes. Note that, similarly 
to our observation on the rectangular driven lattices, the 
latter result implies that dilute swimmer crystals with a 
rectangular symmetry are marginally stable with respect 
to finite amplitude stretch deformations. 



5.2.4 Oblique and face-centered lattices 

The dynamics of active crystals having oblique, or face 
centered symmetries are much more complex. We here 
_briefly highlight some interesting large scale properties, 
:^,nd comment on the stability of these structures. 

The stability matrix of the oblique crystals take a sim- 
ple for the global modes only, q = 0, yet it reveals an 
original dynamics. Indeed for q = we get 



M = 




-ipa - 



ia cos(/3) sin(^) 
os(2^)q . 



(39) 



o 



where f3 is the inclination of the lattice cells and e their 
aspect ratio (Fig. [3]). Recall that the nearest neighbor 
scheme restrains our analysis to weakly tilted and weakly 
anisotropic lattices. Beyond the luq — mode, the other 
two eigenvalues are nonzero, and we obtain 



7r(ea) 



V3/i7[cos(2/3) - cos(6/3)]. (40) 
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For weakly tilted lattices j3 ~ 7r/2, so that the fre- 
quencies are purely imaginary, and thus q = modes are 
unstable. Note that this result does not contradict the sta- 
tionarity of the structure. Indeed, the orientation field is 
here unstable, thereby inducing a coupled translation of 
the lattice, as swimmers rotate. 

Conversely, in the small-g limit, the face-centered lat- 
tices are marginally stable for any amplitude of the hy- 
drodynamic coupling, and phonons and orientation waves 
propagate in a non dispersive manner. For small wave- 
lengths however, and looking specifically at the combi- 
nation (q x = Q,q y — ir/ea), we see that the eigenmodes 
corresponding to compression along the y-direction cou- 
pled to distortions of the orientation grow exponentially 
at a rate ^^xs -y/3p7M- This last results implies that face- 
centered swimmer lattices are unstable for any amplitude 
of the hydrodynamic coupling. 

6 Conclusion 

In this paper we considered theoretically the dynamics and 
stability of both driven and active crystals. With a geom- 
etry of elongated particles under confinement we derived 
the dynamical system quantifying the time-evolution of 
the particle positions and orientations and showed that 
all five planar Bravais lattices are stationary solutions of 
the equations of motion. In the case of particles passively 
driven by an external force we formally demonstrated that 
all five lattices are always marginally stable. The phonons 
modes do not depend on the magnitude of the driving 
force but solely on the orientation and on the symmetries 
of the lattices. We detailed the spatial structure of the 
eigenmodes in the square and hexagonal geometry. 

In the separate case where the particles are actively 
self-propelling we showed that the stability of the parti- 
cle positions and orientations depends not only on the 
symmetry of the crystals but also on the perturbation 
wavelengths and the volume fraction of the crystal. We 
obtained that the square and rectangular lattices are lin- 
early unstable at short wavelengths provided the volume 
fraction of the crystals is high enough. As a difference, 
hexagonal, oblique, and face-centered crystals are always 
unstable. 

The results of our work can be compared with past 
theoretical studies. In the driven case, planar crystalline 
arrangements, were shown to be hydrodynamically unsta- 
ble in a three-dimensional fluid at long wavelengths [27] , 
The results in our paper demonstrate that confinement 
of the crystals, which algebraically screens hydrodynamic 
interactions between the particles, leads to a qualitatively 
different behaviors and all lattices solely support phonon 
modes. 

In the active case, previous work demonstrated the 
presence of long wavelength instabilities in orientation, 
density and stress (see |7l8j and references therein). In 
this past work, aligned suspensions for both pusher and 
puller swimmers were shown to be unstable in the dilute 
regime, and so are isotropic suspensions of pushers [10) 



whereas isotropic puller suspensions, which are linearly 
stable at zero volume fraction, were shown numerically 
to be unstable at high volume fraction |28j . In our pa- 
per, again because of hydrodynamic screening, the stabil- 
ity characteristics of confined active crystals were found 
to be independent of the pusher vs. puller nature of the 
self-propelled particle - the only flow singularity dictating 
hydrodynamic interactions in this case is the potential flow 
dipole whose sign is set by the swimming direction only. 
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